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ABSTRACT 



Context. We detail and benchmark two sophisticated chemical models developed by the Heidelberg and Bordeaux astrochemistry 

groups. 

Aims. The main goal of this study is to elaborate on a few well-described tests for state-of-the-art astrochemical codes covering 

a range of physical conditions and chemical processes, in particular those aimed at constraining current and future interferometric 

observations of protoplanetary disks. 

Methods. We consider three physical models: a cold molecular cloud core, a hot core, and an outer region of a T Tauri disk. Our 

chemical network (for both models) is based on the original gas-phase osu_03_2008 ratefile and includes gas-grain interactions and 

a set of surface reactions for the H-, 0-, C-, S-, and N-bearing molecules. The benchmarking is performed with the increasing 

complexity of the considered processes: (1) the pure gas-phase chemistry, (2) the gas-phase chemistry with accretion and desorption, 

and (3) the full gas-grain model with surface reactions. Using atomic initial abundances with heavily depleted metals and hydrogen 

in its molecular form, the chemical evolution is modeled within 10' years. 

Results. The time-dependent abundances calculated with the two chemical models are essentially the same for all considered physical 

cases and for all species, including the most complex polyatomic ions and organic molecules. This result however required a lot of 

elforts to make all necessary details consistent through the model runs, e.g. definition of the gas particle density, density of grain 

surface sites, the strength and shape of the UV radiation field, etc. 

Conclusions. The reference models and the benchmark setup, along with the two chemical codes and resulting time-dependent 

abundances are made publicly available in the Internet. This will facilitate and ease the development of other astrochemical models, 

and provide to non-specialists a detailed description of the model ingredients and requirements to analyze the cosmic chemistry as 

studied, e.g., by (sub-) millimeter observations of molecular lines. 

Key words, astrochemistry - stars: formation - molecular processes - ISM: molecules, abundances 



^ ■ 1. Introduction 



Astrochemistry plays an important role in our understand- 
ing of the star- and planet-formation processes across 
the Universe (e.g., i Berginet al., 2007'; Kennicutt, '1998'; 
ISolomon & Vande n Bout. 2005; van Dishoeck & Blake, 1998) . 
Molecules serve as unique probes of physical conditions, evolu- 
tionary stage, kinematics and chemical complexity. In astrophys- 
ical objects molecular concentrations are usually hard to pre- 
dict analytically as it involves a multitude of chemical processes 
that almost never reach a chemical steady-state. Consequently, 
to calculate molecular concentrations one has to specify initial 
conditions and abundances, and solve a set of chemical kinetics 
equations. 

Since the first seminal studies on c hemical model- 
ing of the interstellar medium (ISM) by iBates & Spitzej 



1951 



iHerbst & Klempered (Il973h . and IWatson & Salpeteii 
19721), numerous ch e mical models of molecular clouds 
fT992h . 
shells 



(e-g-. 



lAikawa et al 



Hasegawa et al 
200: 



protoplanetary disks 
around late-type stars 



(e.g.. 



Agundez & C ernicharol l2006h . AGN to ri (e.g., |M eiierin k et all, 
i2007 ). and planetary atmospheres (e.g.. [Gladstone et al.L 119961) 
have been developed. These models are mainly based on 
three widely applied ra tefiles: the U DFA (U mist Database 
For Astrochemistr^ d Millar et all Il997t iLe Teuff et all 
2000; Woodalletal., 2007), the OSU database (Ohio State 
Universitjo) ( Smit h et al.L |2004|) . and a network incorporated 
in the "Photo-Diss ociation Region" code from MeudorQ 
(iLe Petit et all l2006h . Recently another database of chemical 
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reactions for the interstellar medium and planetary atmospheres, 
KIDAQ has been opened to the astrochemical community. 
The main aim of KIDA is to group all kinetic data to model 
the chemistry of these environments and to allow physico- 
chemists to upload their data to the database and astrochemists 
to present their models. Several compilation s of surface re- 
action s have als o been presented (e.g., Allen & Robinson, 
1977t JTielens & Hagenl 119821: Hasegawaet al., 1992; 



Hasegawa & Herbstl ll 993t iGarrod & Herbstl l2006h . The 



detailed physics and chemistry processes incorporated in the 
modern theoretical models allow to predict and, at last, to fit 
the observational interferometric data such as, e.g., molecular 
line maps of protoplan etary disks a n d mol e cular clouds (e.g. 



120071: Goicoechea et al.L 



2OO9I: iHenning et al.L 1201 Obh . 



Semenov et al., 2005; iDutrev et al. , 
2009i: iPanic & Hogerheiidei 
However, with the advent of forthcoming high-sensitivity, 
high-resolution instruments like ALMA, ELT, and JWST, the 
degree of complexity of these models will have to be increased 
and their foundations to be critically re-evaluated to match the 
quality of the data (e.g., 'A sensio Ramos et all 120071 ; iLellouchi 
[2OO8; Semenov et al., 2008). 

There are other difficulties with cosmochemical models. 
Apart from often poorly known physical properties of the object, 
its chemical modeling suffer from uncertainties of the reaction 
rates, of which only ^ 10% have been accurately determined 
(see, e.g., Herbst, 1980; Wakelam et al., 2006; Vasvunin et al., 
l2008h . In contrast to these internal uncertainties, another major 
source of ambiguity is the lack of detailed description of chem- 
ical models that often incorporate different astrochemical rate- 
files, initial abundances, dust grain properties, etc., thus mak- 
ing results hard to interpret and compare. It is thus essential to 
perform consistent benchmarking of various advanced chemical 
models under a variety of realistic physical conditions, in par- 
ticular those encountered in protoplanetary disks. To the best 
of our knowledge, the only important benchmarking study at- 
tempted so far focu sed on several PDR physico-chemical codes 
(lRoUigetalil2007l) . 

The ultimate goal of present work is i) to provide a detailed 
description of our time-dependent sophisticated chemical codes 
and ii) to establish a set of reference models covering a wide 
range of physical conditions, from a cold molecular cloud core to 
a hot corino and an outer region of a protoplanetary disk. In com- 
parison with the PDR benchmark, our study is based on a more 
extended set of gas-grain and surface reactions. All the relevant 
software, figures, reaction network, and calculated abundances 
are available in the Internefl 

The organization of the paper is the following. In Section|2l 
the two codes and the chemical network are presented in de- 
tail, along with our approach to calculate the reaction rates. In 
Section [3] we describe various benchmarking models, the phys- 
ical conditions, and the initial chemical abundances. The bench- 
marking runs and their results are presented and compared in 
SectionlH Final conclusions are drawn in Section|5j 



groups. Both codes have been intensively utilized in various 
studies of molecular cloud and protoplanetary disk chemistry, 
e.g. the influence of the reaction rate uncertainties on the re- 
sults of chemical modeling of cores (Vasvunin et al., 2004; 
IWakelam et all l2005l I2OO6) and disks (.Vasvunin et al.- ,20081) . 
modeli ng of the disk chemi c al evolution with tu rbulent dif- 
fusion (Seme nov et al.L 2006; Hersant et al., 7009), interpreta- 
tion of interferometri c data (Semenov et al., 2005; Dutrey et aD , 
120071: ISchrever et all I2O O8; Henning et al., 201()i), and predic- 
tions for ALMA ('Semenov et al.i . i2008). Both codes are op- 
timized to model time-dependent evolution of a large set of 
gas-phase and surface species linked through thousands of gas- 
phase, gas-grain and surface processes. We found that both these 
codes have comparable performance, accuracy, and fast conver- 
gence speed thanks to the use of advanced ODE and sparse ma- 
trix linear system solvers. 

The Heidelberg "ALCHEMIC" code is written in Fortran 77 
and based on the publicly available DVODPK (Differential 
Variable-coefficient Ordinary Differential equation solver with 
the Preconditioned Krylov method GMRES for the solution 
of linear systems) ODE packagq^- The full Jacobi matrix is 
generated automatically from the supplied chemical network 
to be used as a left-hand preconditioner. For astrochemical 
models dominated by reactions involving hydrogen, the Jacobi 
matrix is sparse (having ^ 1% of non-zero elements), so to 
solve the corresponding linearized system of equations a high- 
performance sparse unsymmetric MA28 solver from the Harwell 
Mathematical Software LibrarjQ is used. As ratefiles, both the 
recent OSU06 and osu.2008 and the RATE 06 networks can be 
utilized. A typical run with relative and absolute accuracies of 
10"^ and 10"'^ for the full gas-grain network with surface chem- 
istry (~ 650 species, ~ 7 000 reactions) and 5 Myr of evolution 
takes only a few seconds of CPU time on the Xeon 2.8GHz PC 
(with gfortran 4.3 compiler). 

The Bordeaux "NAUTILUS" code is written in Fortran 90 
and uses the LSODES (Livermore Solver for Ordinary 
Differential Equations with general S parse Jacob ia n ma - 
trix) solver, part of ODEPACK0 toindmarshi Il983h . 
"N AUTILUS" is adapt e d from the original gas-grain model 
of iHasegawa & Herbsa (Il993h and its subsequent evolutions 
made over the years at the Ohio State University. The main 
differences with the OSU model rely in the numerical scheme 
and optimization ("NAUTILUS" is about 20 times faster), the 
possibility to compute ID structures with diffusive transport and 
the adaptation to disk physics and chemistry. The full Jacobian is 
computed from the chemical network without preconditioning. 
For historical reasons, only the OSU rate files are being used, 
although minor adjustments could permit to extend the model 
to other networks. Similar performances with "ALCHEMIC" 
are achieved, the same typical run of the full gas-grain network 
taking a few seconds on a standard desktop computer. 



2. Chemical Models 

2. 1 . Heidelberg and Bordeaux chemical codes 

In this study, we compare two chemical codes, "ALCHEMIC" 
and "NAUTILUS", developed independently over the last sev- 
eral years by the Heidelberg and Bordeaux astrochemistry 



http : //kida . obs . u- bordeaux 1 . f r] 



Below we detail how we model various gas-phase and sur- 
face reactions. 
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2.2. Modeling of Chemical Processes 

Chemical codes solve numerically the equations of chemical ki- 
netics describing the formation and destruction of molecules: 



dn 



~Tr = /_j ^'"'"'"^ ~ "' Zj ^'"' "*" ^^""'"^ ~ ^^ 



dn^ 



Hi 



l,m 



it! 
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(2) 
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where n, and n' are the gas-phase and surface concentrations 
of the j-th species (cm~^), kim and ki are the gas-phase reaction 
rates (in units of s ' for the first-order kinetics and cm^ s ' for 
the second-order kinetics), k^'^'^ and k^'^^ denote the accretion and 
desorption rates (s^^), respectively, and A:''^^ and A;* are surface 
reaction rates (cm^ s~'). 



2.2.1. Gas-phase reactions 

For the benchmarking purposes, we utilize a recent osu_03_2008 
gas-phase ratefilq^- This ratefile incorporates the data for the 456 
atomic, ionic, and molecular species involved in the 4389 gas- 
phase reactions. The corresponding reaction rates are calculated 
as follows, using the standard Arrhenius representation: 



''^^-"[^lA-j\ 



(3) 



where a is the value of the reaction rate at the room temperature 
of 300 K, the parameter p characterizes the temperature depen- 
dence of the rate, and y is the activation barrier (in Kelvin). We 
utilize this expression for all gas-phase two-body processes, e.g. 
ion-molecular and neutral-neutral reactions. 

For the cosmic ray and FUV ionization and dissociation the 
following prescriptions have been used. 



Cosmic ray (CR) ionization & dissociation 



(4) 



where ^cr =13 10 '^ s ' is the adopted CR ionization rate. 
In all benchmarking models any additional sources of ionization 
like the X-ray stellar radiation and the decay of short-living ra- 
dionucleides (e.g., ^''Al and "^''K) are not considered. This small 
set of reactions includes ionization of relevant atoms and dis- 
sociation of molecular hydrogen. The same expression is used 
to compute photodissociation and ionization by the CR-induced 
FUV photons. Also, the CR- and UV-driven dissociation of sur- 
face species is calculated by the Eq.|4](616 processes). 



FUV pliotodissociation & ionization To calculate photoreac- 
tio n rates through the e nvironment, we adopt pre-computed fits 
of lvanDishoeckl ([1988) for a ID plane-parallel slab, using the 
Draine FUV IS radiation field. Unlike in the PDR benchmarking 
study, the self-and mutual-shielding of CO and H2 against UV 
photodissociation are not taken into account for simplicity. The 
corresponding rate is calculated as: 



^Fuv = a exp (-yA^f)x, 



(5) 



where the Ay is the visual extinction (mag.) and the x is the 
unattenuated FUV flux expressed in units of the FUV interstel- 
lar radiation field xa of Draine ( 1978). We do not take the Lyo- 
radiation into account as it requires sophisticated modeling of 
the radiat ion transport - a full - scale benchmarkin g study by it- 
self (e.g.. lPascucci et aUl2004tlRonig et aliboOTh . 



2.3. Gas-grain interactions 

For simplicity, we assume that the dust grains are uniform spher- 
ical particles, with a radius of Og = 0.1 jum, made of amor- 
phous olivine, with density of p^ = 3 g cm"^ and a dust-to- 
gas mass ratio nid/g = 0.01. The surface density of sites is 
Ns = 1.5 10'^ sites cm"^, and the total number of sites per such 
a grain is S - 1.885 10''. The dust and gas temperatures are as- 
sumed to be the same. 

Gas-grain interactions start with the accretion of neutral 
molecules onto dust surfaces with a sticking efficiency of 100% 
(195 processes). Molecules are assumed to only physisorb on 
the grain surface (by van der Waals force) rather than by form- 
ing a chemical bond with the lattice (chemisorption). The rate of 
accretion of a gas-phase species / (cm""* s"') is given by: 



^acc(0 = k.^cc(i)n(i). 



(6) 



where n(i) is the density of gas-phase species / (cm^^), and 
^acc(0 = o-d{v{i))nci is the accretion rate. Here ad - nai 
is the geometrical cross section of the grain with the radius 
ag, Hd is the density of grains (cm"-') and (v(i)) is the ther- 
mal velocity of species / (cms"'). The latter quantity is ex- 
pressed as yj%kBT linjiiCjinp), with T being the gas tempera- 
ture (K), nip - 1.6605410"^"' g is the proton's mass, yu(/) is 
the reduced mass of the molecule / (in atomic mass units), and 
kg - 1.3805410""' ergK"' is the Boltzmann's constant. All 
constants are summarized in Table 2. 

In addition, electrons can stick to neutral grains, producing 
negatively charged grains. Atomic ions radiatively recombine on 
these negatively charged grains, leading to grain neutralization 
(13 reactions in total). The corresponding two-body reaction rate 
is calculated as: 



k{T) 



.300/ n/ 



(7) 



where «h is the total hydrogen nucleus density (cm ^). This 
quantity is calculated by the following expression: 



«H = 



Unip 



(8) 



where p is the gas mass density (g cm"^), and fi - 1 .43 is the 
mean mass per hydrogen nucleus. Consequently, the density of 
grains n^ is expressed as: 



rid 



Al^nalpdli 



(9) 



We consider neither interactions of molecular ions with grains 
nor the photoelectric effect leading to positively charged grains. 
In our simplified model used for benchmarking purposes, the 
surface molecules can leave the grain by only two mechanisms. 
First, they can desorb back into the gas phase when a grain is 



hit by a relativistic iron nucleus and heated for a short while 
|http: //www .physics -Ohio- state.edu/- eric/research .html up to several tens of Kelvin (160 reactions, cosmic -ray induced 



D. Semenov et al.: Benchmarking of chemical models. (RN) 



desorption). Second, in sufficiently warm regions thermal des- 
orption becomes efficient. The thermal desorption for the i-th 
surface molecule is calculated by the Polanyi-Wigner equation; 



kdes(Td) = v(i) exp ( — ^ ). 






- is the characteristic vibrational frequency 



Here v(i) 

of the !-th species, Zsjes is its desorption energy (Kelvin), m the 
mass of the species and Tj is the grain temperature. Desorption 
energies Edes are taken from lGarrod & Herbs J(l2006h . We do not 
distinguish between various thermal evaporation scenarios for 
different molec ules, e.g. via "volcanic" or multilayer desorption 
(ICohings et al.. 2004) . 

The cosmic ray desorption rate is computed as suggested by 
Eq. 15 in'Hasegawa & Herbsll (Il993h . It is based on the assump- 
tion that a cosmic ray particle (usually an iron nucleus) deposits 
on average 0.4 MeV into a dust grain of the adopted radius, im- 
pulsively heating it to a peak temperature T^-p = 70 K (see also 
iLegeret aUll985h . The resulting rate is: 



^crd = /^des(70 K) 



(11) 



kcrd is the thermal evaporation rate at 70 K times the fraction of 
time the grain temperature stays close to 70 K. The fraction / is 
a ratio of a grain cooling timescale via desorption of molecules 
to the timescale of subsequent heating events (e.g., for a O.lyum 
silicate particle and the standard CR ionization rate 1 .3 10 '^ s ' 
these timescales are ~ 10"^ s and ~ 3 10'^ s, respectively). Other 
non-thermal desorption mechanisms that can play an important 
role in chemistry of protoplanetary disks like photodesorption 
are not considered. 



2.4. Surface Reactions 

Surface reactions (532 in the network) are treated in the standard 
rate approximation, assuming only Lan gmuir-Hinshelwood for - 
mation mechanism, as described, e.g., in lHasegawa et al.l(ll992h . 
Surface species are only allowed to move from one surface site to 
another by thermal hopping. When two surface species find each 
other, they can recombine. We assume that the products do not 
leave the surface as the excess of energy produced by such a re- 
action will be immediately absorbed by the grain la ttice, i.e. we 
did no t include the desorption process proposed bv lGarrod et al.l 
(I2007h . 

Rate coefficients for surface reactions between species / and 
j are expressed as follows: 



kj = ^(^diff(0 + Rdm(j))/nd. 



(12) 



Here P is the probability for the reaction to occur. This param- 
eter is 1 for an exothermic reaction without activation energy 
and 1/2 if the two reactants are of the same type. For exother- 
mic reactions with activation energy Eg (or endothermic reac- 
tions), this probability is ffexp ( — j^), with a the branching ratio 
of the reaction (a - 1/3 if there are three reaction channels). 
In some cases, tunneUng effects can increase this probabilit}{3- 
When this is the case, P is co mputed through the following for- 
mula jHasegawa et aU[l992h : 



with b the barrier thickness ( 1 A) and fi the Planck constant times 
27r (1.05459x10"" ergs-'). 

The thermal diffusion rate of species / is: 



(10) ^diff(0 = v(0 exp( — ^^)/5. 



Td 



(14) 



Here ksTdiff is the activation energy of diffusion for the i-th 
molecule. 

The diffusion and desorption energies o f a limited set of 
molecular species have been derived ( e.g., iKatz et all 119991; 
iBisschop et"an.l2006HOberg et al.Ll2007h . Moreover, these ener- 
gies depend on the type of the surface l attice and its structural 
properties (porosity, crystallinity) (e.g., lAcharvva et all 120071; 
Throwe r et al., 2009). While diffusion energy of a species must 
be lower than its desorption energy, the exact ratio is not well 
constraine d for a m ajority of molecule s in our netwo rk (except 
forH2, see lKatz et al.. 1999). FoUowing lRuffle & He rbst (2000), 
we adopted the Tdiff/Tdes ratio of 0.77 for the all relevant species 
in the model. 

It has been proposed that light surface species like H, H2 
and isotopes can also quantum tunnel through a potential well 
of the surface site and thus be able to quickly scan the sur- 
face, greatly increasing the efficiency of surface recomb i nation 
even at very low tem peratures (e.g., iDulev & Williamsl 119861; 
JHasegawa et al.Lfl992h . However, according to the theoretical in- 
terpretation of Katz et al. (1999) tunneling of atomic hydrogen 
on amorphous surfaces does not occur. Therefore, it is not con- 
sidered in our benchmarking study. 

Finally, at very low influx of reacting species having a high 
recombination rate, concentrations on a grain can become very 
low, leading to a stochastic regime. In this case surface chem- 
istry cannot be reliably described by the standard rate equation 
method, which tends to overestim ate the rates. Oth er approaches 
like Monte Carlo techniques (e.g., Vasvunin et al.Ll2009i) . master 
equa tions (e.g.. Green et al., 200 1]), and mod ified rate equations 



(e.g.. lCasellietal.Lll998HGar rod et al.,'2009') should then be uti- 
lized. As it has been shown by Vasvu nin et al. (2009), this only 
happens for rather dilute, warm gas and very small grains, and 
when quantum tunneling of hydrogen is allowed. Thus, the use 
of the standard rate equation approach to model surface chem- 
istry is fully justified in our benchmarking study. 



2.5. Initial abundances & other parameters 

We use an up-to - date set of elemental abundances from 
I Wake lam & Herbsn (l2008l) . The 12 elemental species include H, 
He, N, O, C, S, Si, Na, Mg, Fe, P and CI. Except for hydrogen, 
which is assumed to be entirely locked in molecular form, all 
elements are initially atomic. They are also ionized except for 
He, N and O (see Table [3]l. All heavy elements are heavily de- 
pl eted from the gas phase similar to the "low metals" abundances 
of ' Lee et aP (Il998h . All grains are initially neutral. Since large- 
scale chemical models with an extended set of surface reactions 
rarely reach a chemical steady state, all benchmarking tests run 
over a long evolutionary time span of 10'^ years (with 60 loga- 
rithmic time steps, starting from 1 year). Both chemical codes 
use the same absolute and relative accuracy parameters for the 
solution, 10""*' and 10"^, respectively. 



1/2-1 



(13) 3. Benchmarking models 



P - aexp[-2(b/Ji)(2kBn(i)mpEa) 

'" This process is not included in ALCHEMIC and thus has not been The physical conditions of the five benchmarking ca s es are 
tested in this benchmark. summarized in Table [T] (see also ISnow & McCallL 120061 : 
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iHassel et all l2008l) . These are chosen to represent reaHstic as- 
trophysical object yet to be relatively simple. We decide to focus 
first on physical models representative of less complex astro- 
physical media: a cold core ("TMCl") and a hot corino ("HOT 
CORE"). The "TMCl" model has a temperature of 10 K, a hy- 
drogen nucleus density of 2 x 1 0"* cm"^, and a visual extinction of 
10 mag. The "HOT CORE" model has a temperature of 100 K, 
a hydrogen nucleus density of 2 x 10^ cm"^, and a visual extinc- 
tion of 10 mag. Both models have FUV IS RF;^' - ^Xo of lDraing 
(0978) and ^CR = l.SlO-'^s-i. 

Protoplanetary disks are more complex objects from chem- 
ical point of view. Basically, an outer disk (r S 50 - 100 AU) 
observable with modern radio-interferometers can be divided in 
three distinct parts from top to bottom. The hot and tenuous disk 
surface (usually called "disk atmosphere") is located at above 
4-6 gas scale heights. In this region disk chemistry is similar 
to that of HII and PDR region, with a limited set of primal 
ionization-recombination processes. There ionization is mainly 
governed by the stellar X-ray radiation and stellar and interstellar 
FUV radiation. Closer to the midplane, a partly X-ray and FUV- 
shielded region called the "warm molecular layer" is located (at 
about 1-3 scale heights). This mildly ionized, dense and warm 
layer is a harbor of complex chemistry, where gas-grain inter- 
actions and endothermic reactions are particularly active. There, 
a plethora of more complex species is formed and reside in the 
gas. Below is the highly dense, cold and dark midplane, where 
most molecules are frozen out onto grains, enabling a variety of 
slow surface processes to be active. 

Instead of simulating disk chemistry in the entire disk, we 
decide to pick up a few representative disk regions with highly 
distinct physical conditions, and within the reach of spatial res- 
olution of modern interferometers. We consider three represen- 
tative layers taken at a radius of 98 AU (which requires a sub- 
arcsecond resolution even for nearby objects). Those are desig- 
nated "DISKl" (the disk midplane), "DISK2" (the warm molec- 
ular layer), and "DISK3" (the disk atmosphere). We take physi- 
cal conditions similar to those encountered in the DM Tau disk, 
for which a lot of high-resolution molecular data are available. 
In our simulations, we adopted the 1 + lD steady-state irradiated 
disk model with vertical temperature gradient that represents the 
low-mass Class II protoplanetary d i sk sur rounding the young T 
Tauri star DM Tau (iD'Alessio et alill999l) . The disk has a radius 
of 800 AU, an accretion rate M = 10"^ Mq yr~', a v iscosity pa- 
rameter a - 0.01, and a mass M ^ 0.07 Mq ( iDutrev et al., 1997; 
IPietu etal.L 12007). The central T Tau star has an effective tem- 
perature T, = 4000 K, mass M. = O.SMq, and radius R^ = 2Rq. 

We assumed that the disk is illuminated by the FUV radiation 
from the central star with an intensity ;^f = 410^0 at the distance 
of 100 AU and by the inters tellar UV radiation with intensity xo 
in plane-para l lel geometry dDraind, 1978t Ivan Dishoecki Il988t 
iBrittain et all 120031: iDutrev et all 120071). The visual extinction 
for stellar Ught at a given disk cell is calculated as: 



Av = 



A^H 



mag 



(15) 



where A^e is the vertical column density of hydrogen nuclei be- 
tween the point and the disk atmosphere. Note that, according to 
our definition, the unattenuated stellar FUV intensity for a fixed 
disk radius is the largest in the midplane and gets lower for upper 
disk heights as the distance to the star increases. 

Consequently, the "DISKl" model is located at 6.768 AU 
above the midplane, has a temperature of 1 1 .4 K, a hydrogen 
nucleus density of 5.413 x 10^ cm"^, a visual extinction toward 
the star of 40.35 mag. and an extinction of 37.07 mag. in the 



vertical direction, and the FUV RF intensity of ;\f, = 428.3;^'o- 
The "DISK2" model is located at 29.97 AU above the mid- 
plane, has a temperature of 45.9 K, a hydrogen nucleus density 
of 2.588 X 10^ cm"-', a visual extinction toward the star of 23.23 
mag. and the vertical extinction of 1.939 mag., and the FUV RF 
intensity of ;i', - 393.2;ifo- Finally, the "DISK3" cell is located 
at 45.44 AU above the midplane, has a temperature of 55.2 K, 
a hydrogen nucleus density of 3.669 x 10^ cm"^, a visual ex- 
tinction toward the star of 1.608 mag. and the vertical extinction 
of 0.217 mag., and the FUV RF intensity of ;^ - 353.5;ifo- The 
remaining parameters of all the models are described above. All 
those parameters are summarized in Table 1 . 

4. Benchmarking Results 

Using these 5 representative benchmarking models and our ex- 
tended chemical network, time-dependent chemistry is calcu- 
lated for the entire 10^ years of evolution. A perfect agreement is 
achieved between the Bordeaux and Heidelberg chemical mod- 
els for all considered physical conditions, all species, and all 
time moments. The results for assorted chemical species rep- 
resenting various chemical families and various degree of com- 
plexity are compared in Figs.[T]-[5] The good agreement is in fact 
not surprising when only chemistry is concerned. Both codes 
are based on the rate equations approach to model chemical pro- 
cesses and use well documented and robust procedures to handle 
a multitude of complex physico-chemical processes (e.g., pho- 
todissociation, cosmic ray desorption, photoprocessing of ices, 
surface reactions). 

This agreement required all the constants, reaction rates, and 
parameters of the physical models between the two astrochemi- 
cal codes to match perfectly. There are many important param- 
eters not always mentioned in description of a chemical model, 
and which may hamper a comparison of results. In what fol- 
lows, we discuss major problems that arose during course of our 
benchmarking study. 

The first priority is to check that both codes have exactly 
the same values of fundamental constants expressed in the same 
physical units, and additional useful constants (e.g., year in sec- 
onds, UV albedo of the dust, etc.). Second, the "standard", well- 
defined physical parameters that can be defined as constants in 
the codes (like the cosmic ray ionization rate, the FUV IS, etc.) 
has to be checked. Our next major obstacle was different defini- 
tion of gas particle density: one group uses pure hydrogen nu- 
cleus density while another uses mass density of the gas (with 
all heavy elements included). After we began to use the same 
units for the FUV radiation field and the same conversion factor 
between A^h and Ay our models showed perfect agreement for a 
pure gas-phase chemical network. 

As a next point, we added gas-grain interactions to our mod- 
els. After grain properties (shape, radius, density, surface density 
of sites, porosity, material, etc.) are fixed, apparently the best ap- 
proach is to add reactions between atomic and molecular ions 
with charged grains, grain re-charging, and electron sticking to 
grains. Next, one has to adopt the value of the sticking coef- 
ficient of other gas-phase species. An unexpected issue at that 
stage came from the fact that "ALCHEMIC" used atomic masses 
including isotopes while "NAUTILUS" used atomic masses of 
major isotope only. Next, it is extremely important for disk 
chemical models under comparison to have similar desorption 
and diffusion energies of surface molecules since gas-grain in- 
teractions and build-up of thick complex icy mantles are pow- 
erful in protoplanetary disks. A substantial portion of time to 
reach the perfect agreement between the two models was spent 
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to properly implement desorption mechanisms. We did not con- 
sider UV photodesorption along with thermal and CRP-driven 
desorption because it would require a proper description of the 
UV radiation transport through the disk model. As we relied on 
the same rate equations approach for surface chemistry, the good 
agreement between the full models was reached soon. Here one 
has to be careful with treatment of homogeneous reactions (i.e., 
involving the same species, e.g. H + H — » H2), as their rates by 
statistical arguments are only half of those for heterogeneous re- 
actions. Last, but not least, for easy comparison of results it is 
important to specify in detail the output format of the simulated 
data. 



5. Conclusion & Summary 

We present in this paper the results of several detailed bench- 
marks of the two dedicated chemical codes for the Heidelberg 
and Bordeaux astrochemistry groups. Both codes are used to 
model time-dependent chemical evolution of molecular clouds, 
hot cores and corinos, and protoplanetary disks. The codes use 
the recent osu_03_2008 gas-phase ratefile, supplemented with an 
extended list of gas-grain and surface processes. A detailed de- 
scription of the codes, along with considered chemical processes 
and means to compute their rates are presented. Five representa- 
tive physical models are outlined, and their chemical evolution is 
simulated. The first case ("TMCl") is relevant to the chemistry 
of a dense cloud core, the second case ("Hot Corino") is similar 
to a warm star-forming region, e.g. an inner part of a hot core or 
corino. The last three models correspond to various disk layers 
at the radius of a; 100 AU, representing distinct chemistry. In 
all these benchmarking test runs, despite a large range of con- 
sidered physical parameters like temperature of gas and dust, 
density, and UV intensity, we found perfect agreement between 
the codes. We however had to compare and bring into agree- 
ment step-by-step a large number of parameters and descrip- 
tion of chemical processes in "ALCHEMIC" and "NAUTILUS". 
Moreover, this benchmark study allowed us to correct some mi- 
nor errors and helped us to become fully confident in our chemi- 
cal simulations and predictions. We hope that this study and de- 
tailed report on all components of the codes and models will be 
helpful for other astrochemical models to improve their quality. 
This is an essential step toward development of more sophisti- 
cated chemical models of the ISM and disks, since ALMA, the 
extremely powerful observational facility, will soon become op- 
erational. With this instrument, the quality of molecular interfer- 
ometric maps will drastically improve, allowing observers and 
astrochemists to work together in order to investigate the chem- 
istry of the planetary-forming zones in disks (r < 30 AU) in 
great detail. 
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Fig. 1. Time-dependent abundances as computed with the Heidelberg (crosses) and Bordeaux (solid line) chemical models for 
the "TMCl" case. The X-axes show time in years. The Y-axes are abundances relative to the total amount of hydrogen nuclei. 
The species names are given in each panel. The prefix "g" denotes surface molecules, "CH40" is methanol, and "G-" stands for 
negatively charged grains. 
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Table 1. Physical models 



Model 


T 


n(H+2H2) 


Av 


xl 




[K] 


[cm-3] 


[mag] 




TMCl 


10 


2(4t 


10 


Xo 


HOT CORE 


100 


2(7) 


10 


Xo 


DISKl 


11.4 


5.41(8) 


37.1 


428.3;ro 


DISK2 


45.9 


2.59(7) 


1.94 


393.2;ro 


DISK3 


55.2 


3.67(6) 


0.22 


353.5;ro 



" FUV, 6 < hv < 13.6 eV, (Draine 1978) 
* A(-B) means Ax lO-'* 



Table 2. Constants and fixed parameters 



Name 



Value 



ks 


1.38054 X 10-'" erg R- 


n 


1.05459 X 10-" erg s-' 


foi 


1.3 X 10-"s-l 


Og 


0.1 yum 


Pd 


3 g cm-3 


mjig 


0.01 


Ns 


1.5 X 10'^ sites cm-2 


S 


1.885 X 10' 


ntp 


1.66054 X lO-^"* g 


(^ 


1.43 amu 


T 

^ crp 
f 


70 K 
3 X 10-" 


b 


lA 


Tdiff/Tics 


0.77 



Table 3. Initial abundances 



Species n{X)/nii 



He 


9.00(-2:i^ 


H2 


5.00(-l) 


C+ 


1.2(-4) 


N 


7.6(-5) 





2.56(-4) 


s+ 


8.00(-8) 


Si+ 


8.00(-9) 


Na+ 


2.00(-9) 


Mg+ 


7.00(-9) 


Fe+ 


3.00(-9) 


P+ 


2.00(-10) 


CL+ 


1.00(-9) 



A(-B) means Ax 10 



